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ABSTRACT 


During the course of a study of space shuttle navigation 
capabilities and requi cements, it became necessary to construct a 
digital simulation of the continuous error of the localizer beam of a 
conventional instrument landing system. A discrete mathematical model 
has been developed which is easy to use on a digital computer. This 
model is a system of difference equations driven by a zero-mean 
gaussian random sequence. The model generates an output random sequence 
which is equivalent, for simulation purposes, to the desired random 
process. The equivalence is created by requiring that the first and 
second order statistics of the discrete model be equal to the corre- 
sponding expectations of the random process to be simulated at all 
instants of time which occur in the simulation. 
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DIGITAL SIMULATION OF CONTINUOUS ERROR MODELS WITH APPLICATION 
TO AN INSTRUMENT LANDING SYSTEM ERROR 

Robert B. Merrick and Gerald L. Smith 
Ames Research Center 


SUMMARY 


During the course of a study of space shuttle navigation capabilities 
and requirements , it became necessary to construct a digital simulation of 
the continuous error of the localizer beam of a conventional instrument 
landing system. A discrete mathematical model has been developed which is 
easy to use on a digital computer. This model is a system of difference 
equations driven by a zero-mean gaussian random sequence. The model 
generates an output random sequence which is equivalent, for simulation 
purposes, to the desired random process. The equivalence is created by 
requiring that the first and second order statistics of the discrete 
model be equal to the corresponding expectations of the random process 
to be simulated at all instants of time which occur in the simulation. 

The theory is presented in several distinct sections. Then the 
details of a localizer beam numerical example are given using the same 
format. 


INTRODUCTION 


While participating in a study of space shuttle navigation 
capabilities and requirements, it became necessary to construct a mathe- 
matical model for the digital simulation of the continuous random error 
of the localizer beam of a conventional instrument landing system. The 
localizer beam error random process is nonstationary since the rms error 
decreases with range; however, it has been assumed (ref. 6) that the 
frequency distribution characteristics of this error are described by a 
rational power spectral density function which does not change with time. 
Although a method is well-known (e.g., ref. 1) for obtaining a continuous- 
time representation of a stationary markovian random process from a given 
rational power spectral density function, it is not so obvious how to 
obtain a suitable discrete-time representation for easy digital computer 
implementation. The purpose of this paper is to present a straightforward 
computationally efficient method for constructing such a model and show 
that it may be used for the practical problem at hand. The general 
approach is one suggested by Dr. Stanley F. Schmidt. 
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NOTATION 

COLUMN VECTORS (lower case, underlined) 

Xc an n component differential system response, the state vector 

xj a random sequence which is statistically equivalent to Xc 

at the discrete time points t^ 

y m elements, a linear function of x, y = Hx 

u. r elements--each element is a white noise random sequence 

w n elements--each element is a white noise process 

MATRICES (capital lelers) 

n x n, symmetric 

n x n, a known function of time 

n x n, a known function of time 

m x n, a known function of time 

the expected value of the matrix, x.c xj» the covariance 
matrix of x^. 

the expected value of the matrix, xj 
the expected value of the matrix, w w~ 

the matrix which transitions the state vector from one time 
to another 

n x r, to be determined 
the transpose of the matrix ( ) 


GENERAL 
a 



s 


a constant 

the laplace operator 
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t time 

t a particular time 

AT time interval, t - t 0 or t^ - t(^_ 1 

SUBSCRIPTS 

c continuous 

d discrete 

k an integer 

o i n i I; i .i I 

s steady state 
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THEORY 


1 . Preliminary Considerations 

The motivation for this theoretical development was the need for a 
digital computer algorithm to model the continuous error of an instrument 
landing system localizer beam. It was desired that the algorithm be 
capable of generating a random sequence statistically equivalent to the 
random process to be simulated, and that the algorithm be reasonably simple. 

Unless the given random process happens to be gauss-markov, or some 
other type completely characterized by only a few parameters, the objec- 
tives stated above tend to be mutually contradictory. Therefore, a 
compromise was employed which is common In engineering applications; 
namely, that only first- and second-order statistical equivalence would 
be required. Specifically, at each time point for which the random 
sequence Is defined, Its first- and second-order statistics were required 
to be equal to those of the random process being simulated: 

E<JSd(t k )) = ECXcUk)), EfXdCtkJxJtti)) = Etxcd^xjtt!)} 


The justifications for this compromise are: 

(1) The description of a physical random process if often limited 
to approximating its first- and second-order statistics (such 
Is the case with the example used here), 

(2) First- and second-order statistics are commonly the most 
significant probabilistic attributes of the process, 

(3) The well -developed theory of gauss-markov processes and sequences 
can be applied readily, 

These principles and assumptions are applied In the sequel to develop 
the desired algorithm, the steps in the development being summarized as 
follows : 

(a) From a given rational power spectral density function which 
represents (in the frequency domain) the second-order statistics 
of the error random process, develop a differential equation with 
zero-mean white noise input whose solution would produce elements 
of the desired random process. 

(b) Develop from this differential equation a general relationship 
for the steady-state covariance matrix of the state of this 
continuous system, expressed in terms of the coefficients of 
the differential equation. 
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(c) Develop a general expression for the steady-state covariance 
matrix of the state of a discrete system with a zero-mean white 
noise random sequence input in terms of the system coefficients. 

(d) Equate the covariance expressions for the continuous and discrete 
systems and solve for the coefficients of the required discrete 
system. 

The discrete system (i.e., digital computer algorithm) designed by this 
procedure produces a zero-mean random sequence having the desired second- 
order statistical properties. Usually, as in the example considered herein, 
the process to be simulated is a zero-mean process, or by proper problem 
definition can be made so. However, if not, the desired first-order 
statistical properties can be achieved simply by adding to the generated 
random sequence, a sequence which is the mean of the given error random 
process at the defined time points. 

A certain degree of added generality is obtained automatically in the 
development inasmuch as generally an n-state system must be devised in 
order to model a random process of anything more than rudimentary complexity. 
Since the n-state system has n outputs, and since an infinite variety of 
combinations of these outputs is possible, one can model not only the desired 
random process, but at the same time many other related random processes with 
little additional computation. In particular, it will be noted that a 
variety of non-stationary random processes can be modeled simply by multi- 
plying the system output by an arbitrary time-varying factor. That is, if 
xj is the n-component state vector of the system, and H is an (m x n) 
matrix of time-varying coefficients, then 


2d = H *d 


is a vector random sequence whose m components may have arbitrary 
variances at any given time. 
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2. The Covariance Matrix for a Continuous Error Model (a gauss-markov 
random process) 

A state vector formulation of a linear differential system response is 

Xc = F(t)xc + G(t)w (2) 


where Xc is an n element state vector, w is an n element noise 
vector, and the n * n matrices, F and G, are known functions of time. In 
order that we may apply the useful theory that exists for gauss-markov 
random processes, we now require: 


a) w is zero-mean gaussian white noise 

b) (o) has a zero-mean gaussian distribution 

c) E[w xj(o)] = 0 

Then ^ is a representation of a zero-mean gauss-markov random process 
and 

E{w(t)w T (x)} = Rs(t - t) 

where R is a positive semi -definite n * n matrix and 6 is the Dirac 
delta function. 

For this application, we also require that w be stationary so that 
R is a constant. The associated covariance propagation equation is 
(ref. 2 and 3) 


r X( .U) - FP Xc + P Xc F T + f,W, T (3) 

The solution of the state equation is 

Xc(t) = *(t:t 0 )x c (t 0 ) + j* $>(t:u)G(u)w(u)du (4) 

where <t>(t : t 0 ) is the transition matrix (refs. 4, 5) and the expected 
value of the integral is zero. Equation (4) is a continuous representation 
for the- state vector at any time t in terms of initial conditions at any 
time t 0 , the transition matrix, and the disturbing functions. The 
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transition matrix for a linear constant coefficient differential system 
(F = constant) is a function only of the time difference, t - t 0 ; it is 
not a function of time, t. 

It is readily determined that 

P Xr (t) = «(t:t 0 )P x (t 0 )$ T (t:t 0 ) + B (5) 

V 

where B is a symmetric n x n matrix and 

B = El/ 1 " $(t:u)G(u)w(u)du(j£ $(t:v)G(v)w(v)dv) T } (6) 

t 0 ^ 

It is ordinarily also true (ref. 3) that 

B = <f’(t: v)G(v)RG T (v)$ T (t:v)dv (7) 

to 

Here the integrand is positive semi -definite. Consequently, the integral 
B is also positive semi -definite. 

In any application, the covariance matrix of the continuous random 
process must be explicitly determined. Accordingly, an expression for 
the integral B is needed. The integration is tedious and may be avoided 
if, as in most applications, the random process to be modeled is independent 
of time. When this is true the stable system, which generates the process, 
is in a steady-state condition; the matrices GRGl and F are constant, and 
the transition matrix is a function of the time interval At. In the 
steady-state condition, equations (3) and (5) become: 

0 = FP Xc (t s ) + P Xc (t s )FT + GRGT (8) 

P Xc (t s ) = <KAt)P Xc (t s )$ T (At) + B (9) 

where the notation t s means any time such that the system is in a 
steady-state condition. 

We see that equation (8) determines the elements of the steady-state 
continuous covariance matrix in terms of known constants and equation (9) 
defines B in terms of the P matrix and the transition matrix. Thus, 

B has been determined without integration. 

Continuing the general development and using equations (1), (2), and 
(5) , we write; 
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E(^) = 0 (10) 

Py c (t) = H(t)P Xc (t)HT(t) (11) 

P yc (t) = H(t)*(t:t 0 )P Xc (t 0 )*T H T + HBH t (12) 

where Py c (t) is defined by 

Py c (t) S (13) 


Equations (5), (7), and (12), present the statistical description of the 
random process which is to be simulated with difference equations and 
discrete random variables. 


3. The Covariance Matrix for a Discrete Error Model (a gauss -markov random 
sequence) 

We wish to generate an xd^k) suc ^ that Pxd^k) = Pxc^) at 1 times 
t k , occurring in the digital simulation. We start with the continuous 
equation (4) which has the following discrete analog 


x<j(t k ) = *(t k :t k _ 1 )xd(t k _ 1 ) + Y (t k )u(t k ) (14) 


where Xd 1S an n component vector whose elements are random sequences 
with an initial gaussian distribution. Here u_ is an r component 
vector of gaussian-random sequences with the following characteristics: 


E[uj ( t k ) ] = 0 = E Cu i (t k )x d (t k _ 1 )] 
E[u^(t k )] = 1 for i = 1 to r, r > 1 
E [u-j (t k )uj(t k )] = 0 if i f j 


Consequently, 

E[u(t k )u T (t k )] = I (15) 

The vector ujt k ) is easily obtained from a unit variance gaussian-random 
number generator. 
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Now ^ is a representation of a gauss -markov sequence and 

Px d U k ) ■ *Px d (tk-i)* T + Y(t k )r T (t k ) (16) 

which is the discrete covariance relationship corresponding to the 
continuous relationship expressed in equation (5). 

If we wish to model a random process, which is a linear function of 
a differential equation response as suggested in equation (1), we 
introduce: 

^(t) = HCtJxdCt) (17) 

Py d (t) = H ( t ) P Xd ( t ) hT ( t ) (18) 

and find that 


P yd (t) = H$P Xd (t 0 )$T H T + Hyy T H t (19) 


4. Requirements for the Equivalence of the Continuous and Discrete 
Covariance Matrices 


The two equations which define the manner in which the covariance 
matrices propagate, are repeated here: 


P Xc (t) = $(t:t 0 )P Xc (t 0 )$ T (t:t 0 ) + B (5) 

p x d (tk) = *( t k : - t k-i) p x d ( t .k-i)* T ( t k :t k-i) + y(tk)Y T (tk) (16) 

An inspection of these relationships shows that the discrete variable 
will be statistically equivalent to the continuous quantity if the 
following three conditions are met: 

(a) 4>(tk :t k _ x ) = *(t:t 0 ) 

(b) The discrete variable can be initialized such that 

p x d (tk-i) = p x c (t 0 ) 

(c) A matrix y can be found such that yy T = B 


( 20 ) 



- 10 - 


The example to be presented will illustrate, in detail, how these conditions 
are satisfied in a particular case. 

The other covariance relationships of equation (19) and equation (12) 
are also seen to be statistically equivalent If the preceding three 
conditions are satisified. 


EXAMPLE 


1 . Preliminary Discussion 

A signal whose power spectral density approximates that of the localizer 
noise (ref. 6) can be generated by passing gaussian white noise through a 
filter as indicated below: 


whi te 

S 


k 

y = localizer 

noise 

1 + 3. OS + 2.25S 2 


noise 


It is assumed that initial transients have disappeared so that the 
statistical expectation of the output of the transfer function is not 
time dependent. The gain k is adjusted as a function of range so that 
the variance of the localizer noise representation, y, is appropriate. 

When white noise is passed through a filter of the form 

1 

1 + (2/a ) S + (l/a) z S z 

then when a = 2/3 the derivative of the output has the statistical 
distribution we desire though its magnitude must be adjusted. 


2. An Explicit Presentation of the B Matrix 

A state vector form of the differential equations which represent the 
localizer system is 


x 

x 





0 

w 2 


y = 



X 
I • 

X 


( 21 ) 

( 22 ) 
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where w 2 is a gaussian white noise process of zero-mean and unit variance 
and the system of equations is analogous to equations (1) and (2). In this 
example the vector of noise sources, w, is characterized by 


so that 


E[w(t)w T (i)] 


0 

0 


0 

1 


<$(t - t) 


R = 


0 

0 


0 

1 


(23) 


If we define 




*(AT) 



4> i 2 
$22 



bn 

bi2 



then equation (8) may be applied to this example to get 


P 1 2 = 0 
P 22 = a 3 /4 

Pn = P 22 /a 2 (24) 


Next, applying equation (9) we have 


bn = (a/ 4 )[l - <j>n - ^a 2 ] 
bn = (a/ 4 ) [-$n$2i " $22$i2 a2 ] 
b 22 = (a/ 4 ) [a 2 - $fi " $ 22 a2 ] 


( 25 ) 
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In Appendix A, the transition matrix for this application is determined 
to be 


$11 = (1 + aAt)e" aAt 
$ 12 = Ate" aAt 
$21 = -a 2 Ate" aAt 

$ 22 = (1 - aAt)e _aAt (26) 

and we may find 

bn = (a/4) - (a/4 ) e“ 2aAt [l + 2aAt + 2a 2 (At) 2 ] (27) 

bi 2 = (a/4)e -2aA t[+2a 3 (At) 2 ] (28) 

b 22 = (a/4)a 2 - (a/4)a 2 e~ 2aAt [l - 2aAt +2a 2 (At) 2 ] (29) 


3. Determination and Existence of y 

The components of B have been explicity presented; for a given a 
and At the numerical values of the elements of B can be calculated. 
The next step in obtaining a solution to equation (14) is to determine 
a y such that 


yy T = B (17) 


A lower triangular form for y is sufficient here 



We now have three algebraic equations in three unknowns which are 
readily solvable. If the matrices, y and B, had n rows and columns, 
then there would be n(n + l)/2 equations in the same number of unknowns 
to solve but their solution is straightforward. Appendix B presents the 
details for matrices with three rows and columns. 
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For the present example: 


Yll = bn 1 / 2 

(31) 

Y21 = b 12 /Yl 1 

(32) 

Y 22 = (b22 " b^/bn) 1 / 2 

(33) 


It is clear that the elements of y do exist if 

2 

b 1 1 > 0 and bjib 22 ” b ^2 0 


The matrix B is positive semi -definite, which implies that: 

2 

bn 0, bnb22 - bi2 0 

If bn = o it follows that b 12 is also equal to zero and in this 


case 


y = 


0 

0 




When the matrix B is less than full rank, then the vector ^ may 
be constructed with fewer elements. 


4. Initialization of the Discrete Variable 

All of the major elements necessary for solving equation (14) have now 
been discussed. The random vector u_ is easily obtainable, a transition 
matrix can be calculated, and the matrix y has been presented for the 
example problem. 

The remaining requirement is to determine an initializing vector x^(o). 
This is simple in the given example since the cross correlation term of 
the covariance matrix, Pi 2 » is zero. We wish to construct initial values 
xd(o) and x^o) such that 
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E[Xd(o)] = P u = a/4 
E[*d(o)] = P 22 = a 3 /4 
E[xd(o)xdCo)] = 0 


Take an output from a unit variance gaussian zero-mean random number 
generator (uncorrelated) and multiply this number by P 1 ^ 2 ; this is Xd(o). 

Take a different element from this same gaussian sequence and multiply it 
by P^ 2 ; this is ^( 0 ). These elements form the initial state vector. 


CONCLUDING REMARKS 


In summary, the example random process is modeled for computation by 
programming the following equations 


Xd(t k ) = *(4t)x d (t k _ 1 ) + y(t k )u(t k ) (34) 

y(t k ) = [0 K]x d (t k ) (35) 

where $(At) is given by equation (26), and y(t k ) is obtained from 
equations (27) through (29) and (31) through (33). Initial conditions for 
equations (34) and (35) are 


x^o) = (a/4) 1 / 2 


1 

0 



and u.(t k ) is a 2-component vector whose elements are obtained from a 
gaussian unit-variance random number generator. 

The problem of simulating a statistically defined continuous error on 
a digital computer has been discussed. A technique for creating difference 
equations which fulfil this purpose has been presented. This mathematical 
model is easy to mechanize and it is accurate in the sense that the vector 
of random sequences obtained from it is statistically equivalent, at the 
discrete times, t k , to the vector of random processes being modeled. A 
digital simulation has been constructed for the random process which describes 
the random error of a conventional instrument landing system localizer beam. 
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APPENDIX A 

A DERIVATION OF THE TRANSITION MATRIX FOR THE EXAMPLE APPLICATION 


N(t) 



1 



1 + (2/a )S + (l/a 2 )S 2 



x(t) 


The corresponding equation is 

(S 2 + 2aS + a 2 )L{x(t}> = a 2 L{N(t)} (A1 ) 


When N(t) is identically zero in the equation above, the solution vector 
x(t n ) is related to the same vector at an earlier time, x(t n _j), by a 
multiplying matrix as below 


x <V ' ♦tyViWVi) 


(A2) 


This appendix will explicitly present the transition matrix, 4>, used 
in this application. 

The unforced portion of the solution of equation (A1 ) is 

x ( t) = e" at x n + te" at (x« + ax 0 ) | 

0 u u \ for N = 0 

x(t) = e" at x Q + te- at (-ax 0 -a 2 x 0 )l 


x(t 2 ) 


1 + a(t 2 - ti) 

(t 2 - ti) 

x(t 2 )_ 


_-a 2 ( t 2 - t x ) 

1 - a(t 2 - tj) 




x(ti) 
X { t ! ) 


This is a continuous solution of the homogeneous equation (N = 0). 
It is valid at any time t 2 given the initial conditions at t lt 
Accordingly, the general transition matrix is 
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where Z = 


> 


x(t n ) 


D(1 + aZ) 

DZ 

x(t n -] ) 

x(t n ) 


-Da 2 Z 

D(1 - aZ) 

J (Vi> 



L J 

- 


t n - t n _ ] and D = e" aZ 
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APPENDIX B 

SHOWING THAT THE SYMMETRIC THREE BY THREE MATRIX, B, CAN BE EXPRESSED 
AS yy T WHERE y is LOWER TRIANGULAR 

The relation to be examined here is 


Yll 

0 

o' 


Yll 

Y21 

Y 31 ' 


'bn 

bi2 

b 1 3 

Y21 

Y22 

0 


0 

Y22 

Y 32 

= 

bn 

b22 

b2 3 

Y31 

Y32 

Y33 


0 

0 

Y33 


r 

O" 
»— * 
to 

b2 3 

b 33 


where there are six unknown elements, y-j j . After matrix multiplication, 
the necessary six independent equations may be found by equating corresponding 
matrix elements of the symmetric matricies. 


2 

Yll 

Y11Y21 

Y11Y31 


"bn 

bn 

bi 3 

Y11Y21 

2.2 
Y21 + Y22 

Y21Y31 + Y22Y32 

= 

bn 

b22 

b 23 

YllY31 

(Y21Y31 + Y22Y32 

X 2 2 2 

)Y31 + Y32 + Y33 


bi 3 

b 2 3 

b 33 


and we find that 

Yll = bll 1 / 2 Y 22 = [&22 " Y21] 1 / 2 

Y21 = ^12/Yll Y32 = [^23 ” Y2lY3l]/Y22 

= 1 3/Y 1 1 Y33 = Cb 3 3 - Y31 " Y32] 1 / 2 


Y31 
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